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ABSTRACT 

The need for accurate geometric and radiometric information over large areas has become increasingly important. Laser al- 
timetry is one of the key technologies for obtaining this geometric information. However, there are important application areas 
where the observing platform has its orbit constrained by the other instruments it is carrying, and so the spatial resolution 
that can be recorded by the laser altimeter is limited. In this paper we show how information recorded by one of the other 
instruments commonly carried, a high-resolution imaging camera, can be combined with the laser altimeter measurements to 
give a high resolution estimate both of the surface geometry and its reflectance properties. This estimate has an accuracy 
unavailable from other interpolation methods. We present the results from combining synthetic laser altimeter measurements 
on a coarse grid with images generated from a surface model to re-create the surface model. 

RESUME 

Le besoin deformations geometriques et radiomdtriques precises couvrant de grandes dtendues devient de plus en plus im- 
portant. L'altimetrie laser est une des technologies principals pour obtenir ces informations gdom^triques. Cependant, il 
est des domaines duplication importants ou la plateforme d’observation a son orbite contrainte par les autres instruments 
qu’elle porte, ce qui limite la resolution spatiale qui peut etre enregistr4e par I'altimfetre. Dans cet article nous montrons 
comment reformation enregistr£e par un des autres instruments commun4ment embarqu£s, une camera photographique a 
haute resolution, peut etre combin^e avec les mesures de I’altimfetre laser pour donner une estimation haute resolution £ la 
fois de la geometrie de surface et de ses propriety de reflectivity. Cette evaluation offre une exactitude in^galee par d autres 
methodes d'interpolation. Nous presentons les resultats obtenus en combinant des mesures synthetiques d’altimetre laser sur 
une grille grossiere avec des images produites ^ partir d’un modeie de surface pour reefer le modeie de surface. 

KURZFASSUNG 

Prazise geometrische und radiometrische Informationen uber grosse Areale ist zunehmend von Bedeutung. Die Laser Altimetrie 
ist eine der Schliisseltechnologien zur Gewinnung dieser Daten. Allerdings ist in wichtigen Anwendungsfallen die Laser Altime- 
trie Messung durch weitere Instrumente behindert und daher die raumliche Auflosung eingeschrankt. In dieser Veroffentlichung 
zeigen wir auf wie die von einer hochauflosenden Kamera (einer fast immer installierten Diagnostik) gewonnenen Bilder mit 
den Daten der Laser Altimetrie kombiniert werden konnen urn eine prazise Bestimmung der Oberflachenform und ihrer Reflek- 
tivitatseigenschaften zu ermoglichen. Diese Art der Oberflachenbestimmung erweist sich einer Splineinterpolationen der Laser 
Altimetriedaten uberlegen. Wir zeigen die Ergebnisse der Oberflachenrekonstruktion aus der Kombmation von synthetischen, 
niedrig aufgelosten Laser Altimetriedaten und Bildern. 

1 INTRODUCTION 

The need for accurate geometric information for a variety of 
problems has grown rapidly in the last decades. These needs 
cover a broad field, from monitoring of environmental changes 
such as the deformation rates of glaciers, to the creation of 
3-dimensional digital city models, and the determination of 
the shapes of asteroids and features on planets. The demands 
with respect to the required accuracy are steadily increasing 
(Rees 1990). 

Laser altimetry systems have been able to respond to these 
demands. However, for many applications, only coarse resolu- 
tion sampling is available. This is especially true for planetary 
and small body observations, where the sampling of the sur- 
face is constrained by the orbit of the sensor, and this orbit 
is often determined by the other instruments carried by the 
spacecraft. 

These other instruments usually include a high-resolution op- 


tical imager. These images have been previously used to in- 
fer a surface reconstruction, solving this inverse problem using 
Bayesian probability theory (Smelyanskiy 2000, Morris 2001). 
The accuracy of the reconstruction of the 3-dimensional sur- 
face depends on the geometric information content of the im- 
ages and on additional prior knowledge. Often images from 
mapping orbits do not contain much geometric information 
as the baseline is very small compared to the distance to the 
surface. 

In this paper we show that a dense surface geometry estimate 
can be made by combining the information from a coarse but 
highly accurate grid of height field points from Laser altimetry 
measurements and the limited geometrical information from 
a set of optical images. The resulting surface estimate (both 
geometry and albedo) has a precision unavailable from other 
interpolation methods. At the same time far fewer images 
are needed for a surface reconstruction than without the data 



from the laser altimetry measurements. 

The calculation is a two step process: Using the images and 
a spline interpolation of the laser altimetry data, an approxi- 
mate albedo field of the surface is inferred. This albedo field 
and the spline interpolated surface are the starting points for 
the Bayesian surface reconstruction. The varying accuracy 
of the height field points is taken into account by assigning 
different uncertainty values to the individual points. The un- 
certainty of the laser altimetry measured points is very much 
lower compared with the interpolated values. This approach 
also offers an easy way to combine measurements with differ- 
ent accuracy. The height field points between the grid points 
of the laser altimetry measurements are updated by the ad- 
ditional geometrical constraints of the optical images. We 
present results of the inference of surface models from simu- 
lated height field grids and aerial photographs. The influence 
of different number of images and varying grid resolution is 
shown. 


ing these synthetic observations is straightforward. We also 
assume that the error in these measurements are known. 

3 A BAYESIAN FRAMEWORK 

In this paper the surface geometry is represented by a trian- 
gular mesh and the surface reflectance properties (albedos) 
are associated with the vertices of the triangular mesh. We 
will consider the case of Lambertian surfaces. We will also 
assume that the camera parameters (position and orientation, 
and internal calibration) and the parameters of the lighting 
are known. It is possible to estimate these parameters in a 
similar Bayesian framework, but it is beyond the scope of this 
paper (Morris 2001, Smelyanskiy 2001). 

Thus we represent the surface model by the pair of vectors 
[z p]. The components of these vectors correspond to the 
height and albedo values defined on a regular grid of points 

[zp] = {( 2 i,pi) , i = £ (gx+py)} q,p = 0, 1, . . . (1) 


2 THEORY 

The objective here is to infer a surface model using the avail- 
able data, in this case, laser altimeter measurements and op- 
tical images. Bayesian inference has, for some time now, been 
the method of choice for many inference problems, enabling 
accurate estimation of parameters of interest from noisy and 
incomplete data (Bernardo 1994). It also provides a con- 
sistent framework for the incorporation of multiple, distinct, 
data sets into the inference process. The general approach is 
illustrated in figure 1. The figure shows that synthetic obser- 
vations of the model are made using a computer simulation 
of the observation process, and that these are compared with 
the actual observations. The error between the actual and 
the simulated observations is used to adjust the parameters 
of the model, to minimize the errors. Bayes theorem tells us 
directly how much weight to assign to the two sources of er- 
rors, those coming from the image measurements and those 
coming from the laser altimeter measurements. 

The surface model we use here is a triangulated mesh. At 
each vertex of the mesh we store the height and the albedo. 
As discussed above, to be able to infer the surface heights 
and albedos, we must first be able to simulate the data that 
would be recorded from the surface. 

Generating images from the surface model is the area of com- 
puter graphics known as rendering (Foley 1990). It is impor- 
tant to note, however, that much recent work in computer 
graphics is unsuitable for our purpose, as it works in image 
space , where the fundamental unit is the image pixel, and any 
given pixel is coloured by light from one and only one surface 
element. This results in artefacts due to the relative sizes of 
the projections of the surface elements onto the image plane 
and their discretization into pixels. These artefacts are par- 
ticularly noticeable along the edges of the surface elements 
(aliasing). For this work we require the renderer to operate 
in object space, and below we will briefly describe such a sys- 
tem. We also note that an object-space renderer can also 
compute the derivatives of the pixel values with respect to 
the surface model parameters. This is crucial in enabling ef- 
ficient estimation of the surface model parameters, and will 
be described in more detail below. 

We are also required to produce synthetic laser altimeter mea- 
surements. We make the approximation that the laser altime- 
ter makes point measurements of the surface, and so produc- 


where i is the elementary grid length, x, y are an orthonormal 
pair of unit vectors in the (x,y) plane and i indexes the posi- 
tion in the grid. The pair of vectors of heights and albedos 
represents a full vector for the surface model 

u = [z p). (2) 


To estimate the values of z,p from the laser altimeter and 
image data, we apply Bayes theorem which gives 

p(z, p |L, 1 1 . . . If ) ocp(L,/i .. . I f\z>p) p(z, P )> (3) 

where L is the laser altimeter data = the 

image data. This states that the posterior distribution of the 
heights and the albedos is proportional to the likelihood - 
the probability of observing the data given the heights and 
albedos - multiplied by the prior distribution on the heights 
and albedos. 

Given the surface description, the images and the laser altime- 
ter measurements are conditionally independent, and equa- 
tion 3 can be written as 


p(z,p\L,h ...If) oc p(L\z,p)p(h ■■ . I F \z, p) p(z, p), 


where we now hsve two independent likelihood terms, one for 
each data stream. 

The prior distribution is assumed to be Gaussian 


p( 2 ,p)ocexp(-luS ‘u r ), 
] ’ 


(4) 


Q/°l 
0 Q/< 


where the vector of the surface model parameters u is de- 
fined in (2). The inverse covariance matrix is constructed to 
enforce a smoothing constraint on local variations of heights 
and albedos. We penalize the integral over the surface of the 
curvature factor c(x, y) = + z\ y + 2 z% yt and similarly 

for albedos. The two hyperparameters ok and a p in equa- 
tion (4) control the expected values of the surface- averaged 
curvatures for heights and albedos. 

This prior is placed directly over the height variables, z, but 
albedos are only defined over the range [0— 1]. To avoid this, 




Figure 1: Outline of the Bayesian approach to surface reconstruction from images and laser altimeter measurements 


we use transformed albedos p[ in the Gaussian (4), where p- t 
are defined by: 

Pi — “ Pi)) » u-+[2/?]. (5) 


In the vector of model parameters u values of p are replaced 
by values of p‘ . 

For both the likelihoods we make the usual assumption that 
the differences between the observed data and the data syn- 
thesized from the model have a zero mean, Gaussian distri- 
bution. So for the laser altimeter measurements we have 


p(L\z,p) oc exp 


2 af 


( 6 ) 


where the summation is over the individual measurement 
points Li, and U{z y p) denotes the laser altimeter measure- 
ments synthesized from the model. The parameter a f is the 
variance of the laser altimeter measurement system. 

We also assume that the images If comprising the data are 
conditionally independent, giving 


p(h ■■ -If\z,p) oc exp 




2a; 


where Ifp(z y p) denotes the pixel intensities in the image / 
synthesized from the model, a\ is the noise variance and the 
summation is over the pixels (p) and over all images (/) used 
for the inference. 

Consider the negative log-posterior. 


C(z,p) oc 


^ f T {Ifv Ifp( z iP)) 

al 

-Li(z } p)) 2 
+ ^2 


( 7 ) 


where x = u - u 0 is a deviation from a current estimate u 0 . 
£ is a nonlinear function of z y p and the MAP estimate is 
that value of z : p which minimizes £(z,p). 

The crux of the problem is thus how to minimize C. We apply 
a gradient method, using an initialization based on a spline 
interpolation of the laser altimeter measurements. 

Making the assumption that the laser altimeter makes point 
measurements of one of the vertices of the mesh, equation 6 
can be written as 

p(L\z,p) oc exp (-1/2(1 - lo)£f x (l - lo) T ) 


where l 0 is the vector of actual laser altimeter observations, 
and I are the corresponding entries taken from the z vector. 
The inverse covariance matrix a^ 1 is a diagonal matrix with 
l/af on the leading diagonal. 

The term for the image measurements is more complex, as 
I(z,p) is the rendering process. To make progress with min- 
imizing C(z,p) we linearize I(z,p) about an initial estimate, 

zo y po 

/(z,p) = /(ib,Po) + Dx, D = j - 7 ^ , j ( 8 ) 


where D is the matrix of derivatives evaluated at 20 , po- Then 
the minimization of L(z,p) is replaced by minimization of the 
quadratic form: 


C! = 
A = 
b = 


-xix-bx, 


S' 1 + 


DD 


(J -/(zo,po)) 
^2 


X = 


D 


u - u 0 , 


( 9 ) 


where tf 1 is now a large square matrix (of dimension 
length(f) + length(p)), where the diagonal elements corre- 
sponding to the vertices for which there are laser altimeter 
measurements take values 1 /af and all other entries are zero. 
The entries of uo corresponding to the laser altimeter mea- 
surements are set to the observed values, and the remaining 
height values are initialized using a spline interpolation. This 
interpolation and the albedo initialization will be described 
below. 

In equation 9 A is the Hessian matrix of the quadratic form 
and vector b is the gradient of the likelihood £ computed at 
the current estimate. We search for the minimum in x using 
a conjugate-gradient method (Press 1992). 

Thus the most difficult part of finding the MAP is the re- 
quirement to render the image and compute the derivatives 
for any values of the surface model parameters. We discuss 
this computation in some detail in the next sections. Here 
it is sufficient to note that while forming / using only object 
space computation (see section 4) is computationally expen- 
sive, we can compute D at the same time for little additional 
computation. Also the derivative matrix is sparse with the 
number of nonzero entries a few times the number of model 
parameters. This makes the process described above a prac- 
tical one. 

The log-posterior is potentially multi-modal, and so it is im- 
portant to begin the optimization from a good initialization. 
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Figure 2: Geometry of the triangular facet, illumination direc- 
tion and viewing direction. z s is the vector to the illumination 
source; z v is the viewing direction. 

In order to do this, the high resolution surface estimation 
proceeds as follows. 

1. Use a spline interpolation of the laser altimeter mea- 
surements to produce an initial height field estimate at 
the desired high resolution. 

Also produce the matrix Ef 1 , with zeros everywhere 
except those diagonal entries corresponding to the 
points of the high resolution surface grid that are mea- 
sured by the laser altimeter. These values are \/o\. 
Initialise all the albedos to 0.5. 

2. Render the surface generated in step 1 and compute 
the derivative matrices D (one for each image). Set 
to zero all the derivatives with respect to the surface 
heights. This fixes the heights at their current values 
in the optimization step (below), so that only the sur- 
face albedo values are inferred. 

Starting from the surface from step 1, and using the 
derivative matrices calculated above, use the conjugate 
gradient algorithm to minimize the linearization of the 
log-posterior in equation 9. This produces a good ini- 
tialization for the final optimization. 

3. Render the surface generated in step 2 and compute 
the derivative matrices D. 

Starting from the surface in step 2, use the conjugate 
gradient algorithm to minimize the linearization of the 
log-posterior. 

4. Repeat step 3 until convergence. 

The minimum found is the final surface estimate, which 
combines the information from the laser altimeter mea- 
surements and the visible images. 

4 FORMATION OF THE IMAGE AND THE 
DERIVATIVE MATRIX. 

The task of forming an image, /, given a surface descrip- 
tion, z,p, and camera and illumination parameters is the 
area of computer graphics known as rendering (Foley 1990). 
Most current rendering technology is focused on producing 
images which are visually appealing, and producing them very 
quickly. As discussed in the introduction, this results in the 


use of image-space algorithms, with the fundamental assump- 
tion that each triangle making up the surface, when projected 
onto the image plane, is much larger than a pixel. This makes 
reasonable the assumption that any given pixel receives light 
from only one triangle, but does produce images with artifacts 
at the triangle edges. Standard rendering also produces inac- 
curate images if the triangles project into areas much smaller 
than a pixel on the image plane, as the pixel will then be 
colored with a value coming from just one of the triangles. 

Clearly this approach is not suitable for high-resolution 3D 
surface reconstruction from multiple images. The triangles 
in a high-resolution surface may project onto an area much 
smaller than a single pixel in the image plane (sub-pixel res- 
olution). Therefore, as discussed in the introduction, for 
our system we implemented a renderer for triangular meshes 
which performs all computation in object space. At present 
we neglect the blurring effect due to diffraction and due to 
the role of pixel boundaries in the CCD array. Then the light 
from a triangle as it is projected into a pixel contributes to 
the brightness of the pixel with a weight factor proportional 
to the fraction of the area of the triangle which projects into 
that pixel. This produces anti-aliased images and allows an 
image of any resolution to be produced from a mesh of ar- 
bitrary density, as required when the system performing the 
surface inference may have no control over the image data 
gathering. 

Our renderer computes brightness I p of a pixel p in the image 
as a sum of contributions from individual surface triangles 
A whose projections into the image plane overlap, at least 
partially, with the pixel p. 

= (10) 
A 

Here is a radiation flux reflected from the triangular facet 
A and received by the camera, and is the fraction of the 
flux that falls onto a given pixel p in the image plane. In the 
case of Lambertian surfaces and a single spectral band $a is 
given by the expression 

$&= p E(a s ) cos a v cos K 6 Afl> (11) 
E(a 3 ) = A (l 3 cosa s +1*) . 

AO = S/d 2 . 

Here p is an average albedo of the triangular facet. Orienta- 
tion angles a* and a v are defined in figure 2. E(a s ) is the 
total radiation flux incident on the triangular facet with area 
A. This flux is modeled as a sum of two terms. The first 
term corresponds to direct radiation with intensity X s from 
the light source at infinity (commonly the sun). The second 
term corresponds to ambient light with intensity X a . The pa- 
rameter 0 in equation. (11) is the angle between the camera 
axis and the viewing direction (the vector from the surface to 
the camera); The exponent k is the lens fa Hoff factor. AH 
in (11) is the spatial angle subtended by the camera which is 
determined by the area of the lens S and the distance d from 
the centroid of the triangular facet to the camera. 

We identify the triangular facet A by the set of 3 indices 
(io 5 ii,i 2 ) from the vector of heights (1) that determines the 
vertices of the triangle in a counterclockwise direction (see 
figure 2). In the r.h.s of equation (11) we have omitted 
for brevity those indices from all the quantities associated 
with individual triangles. The average value of albedo for the 
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Figure 3: The intersection of the projection of a triangular 
surface element (io, ii, 12 ) onto the pixel plane with the pixel 
boundaries. Bold lines corresponds to the edges of the poly- 
gon resulting from the intersection. Dashed lines correspond 
to the new positions of the triangle edges when point Pi 0 is 
displaced by £P 


triangle in (11) is computed based on the components of the 
albedo vector p corresponding to the triangle indices 

Pa = pio>i »>3 “ 3(^0 P* 1 P*a)- (1^) 

We note that using average albedo (12) in the expression for 
is an approximation which is justified when the albedo 
values vary smoothly between the neighboring vertices of a 
grid. 

The area A of the triangle and the orientation angles in (11) 
can be calculated in terms of the vertices of the triangle Pj 
(see figure 2) as follows: 

n • z s — cos a 3 , n • z v == cos a v , (13) 

~ __ Vi Q.»l X V »l>»3 v ■ r P- - P 1 

2A ' ,J J 

Here n is a unit normal to the triangular facet and vectors of 
the edges of the triangle v S j are shown in figure 2. 

We use a standard pinhole camera model with no distortion in 
which coordinates of a 3D world point P = {x : y^z) are first 
rotated with the rotation matrix R and then translated by the 
vector T into camera coordinates, yielding P c = (x C) yc,z c ) 

P C = RP + T (14) 


(R and T are expressed in terms of the camera registration 
parameters (Hartley 2000). We do not give them explicitly 
here). After the 3D transformation given in (14), point P c 
in the camera coordinate system is transformed using a per- 
spective projection into the 2D image point P = ( x,y ) using 
a focal length / and aspect ratio a. 
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Figure 4: The initial synthetic surface model (Duckwater, 
NV). (Note the exaggerated vertical scale.) 


Here A a is the area of the projected triangle on the image 
plane and A po | ygon is the area of the polygon resulting from 
the intersection of the projected triangle and the pixel p (see 
figure 3). 

4.1 Computation of the derivative matrix. 

The inference of the surface model parameters depends on 
the ability to compute the derivatives of the modeled obser- 
vations / with respect to the model parameters. According 
to equation (10), the intensity I p of a pixel p depends on the 
subset of the surface parameters, (heights and albedos), that 
are associated with the triangles whose projections overlap 
the pixel area. 

The derivatives Ip with respect to logarithmically transformed 
albedo values are easily derived from equations (5), (10) and 
( 11 ). 

In our object-space renderer, which is based on pixel-triangle 
geometrical intersection in the image plane, the pixel inten- 
sity derivatives with respect to the surface heights have two 
distinct contributions 


di p 

dz\ 



d$A 

dz\ 


+ 


dfk 

dz\ 


(17) 


Variation of the surface height z\ gives rise to variations in 
the normals of the triangles associated with this height (in a 
general triangular mesh, on average 6 triangles are associated 
with each height) and this produces the derivatives of the to- 
tal radiation flux $a to the camera from those triangles. This 
is the first term in equation (17). Also, height variation gives 
rise to the displacement of the corresponding point which is 
the projection of this vertex on the image plane. This results 
in changes to the areas of the triangles and polygons with 
edges containing this point (see figure 3). This produces the 
derivatives of the fractions /£, the second term in equation 
17. Details of these derivatives can be found in (Smelyanskiy 
2000, Morris 2001, Smelyanskiy 2001). 


We use 2D image projections of the triangular vertices Pj to 
compute the area fraction factors /£ for the surface triangles 
(cf. Eq. (10)) 

fl = (16) 


5 RESULTS 

Figure 4 shows the synthetic surface that we will use to 
demonstrate our methodology. The topography is taken from 
the USGS DEM of Duckwater, Nevada. A LANDSAT-TM 




camera 

(75,150,2000) 

image 1 

look at 

(150,150,0) 


view up 

(0,1,0) 


camera 

(225,150,2000) 

image 2 

look at 

(150,150,0) 


view up 

(0,1,0) 


Table 1: Camera parameters used to generate the images in 
figure 5 


Figure 5: Images of the synthetic surface 

image was co-registered with the DEM, and the values of 
one band were used in place of the true albedos. This results 
in the surface shown. One unit is approximately 180 meters. 

Figure 5 shows two images rendered from the surface, and 
table 1 gives the positions and orientations of the (synthetic) 
cameras. The cameras were positioned to approximate satel- 
lite observations. The two images look very realistic. Note 
that the images appear very similar due to the proximity of 
the two camera positions. There is limited geometric infor- 
mation available from the images alone. 

Figure 6 shows a surface from a grid of 9 x 9 points extracted 
from the surface in figure 4. The major terrain features have 
alt been sampled, but clearly it is a very poor representation of 
the surface. This is taken as the laser altimeter observations 
of the surface. 

Using the images and the 9 x 9 grid, we will now go through 
the surface estimation procedure that was detailed above. 

Figure 7 shows the result of using the standard spline inter- 
polation to expand the 9 x 9 grid to the full resolution of the 
surface. The result is a smooth surface showing the major 
features, but note that it contains no more information than 
the coarse surface. 

Keeping the heights fixed at this surface, we then use the im- 
ages to infer initial values for the albedos. The result is shown 
in figure 8. Note that this is not simply the back-projection 
of the of the images onto the surface, The information in 
both images has been optimally combined to give the albedo 
estimates. This surface is now a passable approximation to 
the original surface, as it has high resolution albedo infor- 
mation providing rich visual detail, but clearly it contains no 
topographic detail. 

Figure 9 shows the final inferred surface. This is much im- 
proved over the surface in figure 8. It shows that much of 
the detail of the topography has been extracted from the data 
and incorporated into the model. The error surfaces show in 
figures 10 and 11 show clearly the improvement in the sur- 
face estimate. These error surfaces show both the height 
error and the albedo error as a shaded surface - the topog- 
raphy shows the height error, and the colour of the surface 



Figure 6: Grid from the 9 x 9 simulated laser altimeter mea- 
surements. 



Figure 7: Spline interpolation of the synthetic laser altimeter 
measurements. 

shows the albedo error. The rms errors for the interpolated 
surface are 0.64 (per vertex) for the heights and 2.6 x 10 
for the albedos. For the final infered surface the errors are 
0.01 for the heights and 1.3 x 10“ 5 for the albedos. Note 
that the albedo values from the initialization are already quite 
good (as can be seen on figure 8), however the inference pro- 
cess produces a topography which is very significantly more 
accurate. 

6 CONCLUSIONS AND FUTURE EXTENSIONS 

We have presented the theory and practice of using Bayesian 
methodology to combine the information in laser altimeter 
measurements and visible images into a single, high resolution 
surface model. We have shown on synthetic data that the two 
data sets can be combined into a single high resolution model 
that is more detailed than could be provided by either data 
stream alone. 

Current work is proceeding towards applying the demonstra- 
tion system to real data, including NASA mission data. Work 
in this area is devoted to sensor modeling (producing the syn- 
thetic images and derivative matrices for the actual imaging 
sensor, rather than an idealization of it), estimation of the 
camera positions to sub-pixel accuracy, better control of the 
smoothness prior on the surface, and better initialization of 
the optimization procedure. 







Figure 8: Interpolated surface with inferred albedos. 



Figure 9: Final inferred surface using the Bayesian approach 
to combining the laser altimeter and visible images 



Figure 11: Error surface for the inferred surface 
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Figure 10: Error surface for the interpolated surface 


